First we load up the packages we’ll be working with
# remotes::install_github("mathesong/bloodstream")
# remotes::install_github("mathesong/kinfitr")
library(tidyverse)
library(bloodstream)
library(kinfitr)
library(job)
theme_set(theme_light())
Here is an optional chunk to download the minified dataset. However, it should already be included in the binder.
download.file("https://www.dropbox.com/scl/fi/ax1t3lfop3pp3h4cir0m7/ds004869.zip?rlkey=6z7e9wopud9ngu26ekbd3kfks&dl=1",
destfile = "ds004869.zip")
unzip(zipfile = "ds004869.zip")
Here we can run a basic blood analysis. For anyone who wants to run it, I’ve placed it in a job call, which runs it in the background, but it may or may not make your session slow.
job({
bloodstream("ds004869/")
})
Alternatively, if you want to use a configuration file, this is how you could do so. You can inspect the config file in the relevant directory. But best not to run this: this takes a little longer, and it is not in a job call either, so you won’t be able to do anything until it is completed.
bloodstream("ds004869/", configpath = "ds004869/code/config_tutorial.json")
We want the arterial input function data, which are called “input”.
bloodstream_data <- bloodstream_import_inputfunctions("ds004869/derivatives/bloodstreamtutorial/") %>%
select(-measurement)
Let’s plot a few of them to see. These contain the interpolated
curves, and can be used with kinfitr to perform kinetic
modelling.
head(bloodstream_data, n = 6)
bloodstream_data %>%
slice(1:6) %>%
pull(input) %>%
map(plot)
[[1]]
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).
[[2]]
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).
[[3]]
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).
[[4]]
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).
[[5]]
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).
[[6]]
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).
Now we want to load up the petprep outputs. There will
eventually be a function for automating the next few steps, but the PET
BIDS Derivatives guidelines are not yet nailed down, so we’ll do it
manually for now.
First we parse the files and load in the data for the tacs and the morphology.
petprep_data <- kinfitr::bids_parse_files("ds004869/derivatives/petprep_extract_tacs/") %>%
unnest(filedata) %>%
filter(str_detect(path_absolute, "gtmseg"))
tacs <- petprep_data %>%
filter(measurement=="tacs") %>%
filter(is.na(pvc)) %>%
mutate(tacdata = map(path_absolute, ~read_delim(.x, delim="\t",
show_col_types = FALSE)))
morphdata <- petprep_data %>%
filter(measurement=="morph") %>%
mutate(morphdata = map(path_absolute, ~read_delim(.x, delim="\t",
show_col_types = FALSE)))
Now we have lots of little regions of interest, but they are all very small. Typically, we would want to combine regions into larger regions that we are interested in. We do so using volume-weighted averaging, which requires both the morphology and the TAC information. This will also soon be automated once the PET BIDS Derivatives are more finalised.
We’ll create a frontal cortex region, a striatal region and a hippocampus-amygdala region.
First, let’s define the regions by their constitutent parts.
frontal_regions <- morphdata$morphdata[[1]] %>%
filter(str_detect(name, "frontal")) %>%
pull(name)
frontal_regions
[1] "ctx-lh-caudalmiddlefrontal" "ctx-lh-lateralorbitofrontal" "ctx-lh-medialorbitofrontal" "ctx-lh-rostralmiddlefrontal"
[5] "ctx-lh-superiorfrontal" "ctx-lh-frontalpole" "ctx-rh-caudalmiddlefrontal" "ctx-rh-lateralorbitofrontal"
[9] "ctx-rh-medialorbitofrontal" "ctx-rh-rostralmiddlefrontal" "ctx-rh-superiorfrontal" "ctx-rh-frontalpole"
striatal_regions <- morphdata$morphdata[[1]] %>%
filter(str_detect(name, "Putamen|Accumbens|Caudate")) %>%
pull(name)
striatal_regions
[1] "Left-Caudate" "Left-Putamen" "Left-Accumbens-area" "Right-Caudate" "Right-Putamen" "Right-Accumbens-area"
hipamg_regions <- morphdata$morphdata[[1]] %>%
filter(str_detect(name, "Hippocampus|Amygdala")) %>%
pull(name)
hipamg_regions
[1] "Left-Hippocampus" "Left-Amygdala" "Right-Hippocampus" "Right-Amygdala"
Now we combine the TACs and the region sizes side-by-side
selected_tacs <- select(tacs, c(ses:rec, tacdata)) %>%
inner_join(select(morphdata, c(ses:rec, morphdata))) %>%
group_by(sub, ses) %>%
mutate(tacdata = map(tacdata, ~pivot_longer(.x,
cols = `Left-Cerebral-White-Matter`:`ctx-rh-insula`,
names_to = "name", values_to = "TAC"))) %>%
mutate(tacdata = map2(tacdata, morphdata, ~inner_join(.x, .y, by="name")))
Joining with `by = join_by(ses, sub, task, trc, acq, run, rec)`
Then we perform the volume-weighted averaging itself. I define a function to do it for one region, and another function to do it for multiple regions.
volume_weighted_average_tac <- function(tacdata, regions, regionname) {
tacdata_combined <- tacdata %>%
filter(name %in% regions) %>%
group_by(frame_start, frame_end) %>%
mutate(volume_tot = sum(`volume-mm3`),
volume_frac = `volume-mm3` / volume_tot,
TAC_frac = TAC * volume_frac) %>%
summarise(!!regionname := sum(TAC_frac),
.groups = "keep") %>%
ungroup()
return(tacdata_combined)
}
volume_weighted_average_tacs <- function(tacdata, regions_list) {
regionnames <- names(regions_list)
out <- map2(regions_list, regionnames, ~volume_weighted_average_tac(tacdata, .x, .y))
out <- purrr::reduce(out, dplyr::inner_join, by = c("frame_start",
"frame_end"))
return(out)
}
regions_list <- list(
Frontal = frontal_regions,
Striatum = striatal_regions,
HippAmg = hipamg_regions
)
selected_tacs <- selected_tacs %>%
group_by(sub, ses) %>%
mutate(selected_tacdata = map(tacdata, ~volume_weighted_average_tacs(.x,
regions_list)))
Now we combine the blood data and PET data so that we can perform modelling using both of them.
modeldata <- selected_tacs %>%
select(ses:rec, tacdata = selected_tacdata) %>%
inner_join(bloodstream_data)
Joining with `by = join_by(ses, sub, task, trc, acq, run, rec)`
Many of these steps will be automated when the PET BIDS Derivatives are finalised. But for now, we’ll do these ourselves.
The TAC data are in Bq/mL, and the bloodstream data are in kBq/mL.
This is because Bq/mL is hard to read, but it is the SI unit. Future
functions within kinfitr will automatically convert
them.
modeldata <- modeldata %>%
mutate(tacdata = map(tacdata, ~.x %>%
mutate(across(.cols = Frontal:HippAmg,
~unit_convert(.x,
from_units = "Bq",
to_units = "kBq")))))
We only have frame_start and frame_end, but we want the midtimes and durations for modelling. We also want the times in minutes for modelling, so we’ll convert them.
modeldata <- modeldata %>%
mutate(tacdata = map(tacdata, ~.x %>%
mutate(frame_start = frame_start / 60,
frame_end = frame_end / 60) %>%
mutate(frame_dur = frame_end - frame_start,
frame_mid = frame_start + 0.5*frame_dur)))
There are many different weighting functions, many of which are
included in kinfitr. I’ll use the default one for now. We
usually run the weights on a large global grey matter brain region. For
today, we’ll just take a mean TAC of our three regions and use that.
modeldata <- modeldata %>%
mutate(tacdata = map(tacdata, ~.x %>%
mutate(meanTAC = rowMeans( .x %>% select(Frontal:HippAmg) )) %>%
mutate(weights = weights_create(t_start = frame_start,
t_end = frame_end,
tac = meanTAC))))
Now we will fit the two-tissue compartment model to a single TAC just to see how it looks.
First we should fi
fit_2tc <- twotcm(
t_tac = modeldata$tacdata[[1]]$frame_mid,
tac = modeldata$tacdata[[1]]$Frontal,
input = modeldata$input[[1]],
weights = modeldata$tacdata[[1]]$weights,
vB = 0.05, multstart_iter = 5)
plot(fit_2tc)
fit_2tc$par
That worked!
Notice that in the parameters, we have estimated an
inpshift. This is the delay: a shift in the timing of the
blood (input) data to the same time as the PET. This is necessary both
because the tracer arrives at the brain at a slightly different time
from when it arrives at the arm where we measure it, but also because it
is hard for these two clocks to be perfectly synchronised.
We would usually fit the delay using the first 5 or 10 minutes of the acquisition using a large region for each individual with a one-tissue compartment model and check that. Then we would use that value when fitting the later models. That information is covered more fully in the full on my website.
Here we’ll use a linearised model because they fit more quickly, in this case Logan. They are also less sensitive to the delay time, so we can probably reasonably safely ignore it. But most linearised models require a t* time to operate.
Let’s choose an appropriate t* value
modeldata %>%
ungroup() %>%
filter(ses=="baseline") %>%
slice(1:5) %>%
mutate(tstarplot = map2(tacdata, input,
~Logan_tstar(
t_tac = .x$frame_mid,
lowroi = .x$HippAmg,
medroi = .x$Frontal,
highroi = .x$Striatum,
input = .y,
vB = 0.05)
)) %>%
pull(tstarplot)
Warning: There were 15 warnings in `mutate()`.
The first warning was:
ℹ In argument: `tstarplot = map2(...)`.
Caused by warning:
! Removed 1 row containing missing values or values outside the scale range (`geom_point()`).
ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 14 remaining warnings.
[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
Ok, let’s use 13 frames, corresponding with 57.5 minutes.
Let’s focus on the frontal cortex, and run the fitting for all measurements.
modeldata <- modeldata %>%
group_by(sub, ses) %>%
mutate(Logan_fit = map2(tacdata, input, ~Loganplot(t_tac = .x$frame_mid,
tac = .x$Frontal,
input= .y,
tstarIncludedFrames = 13,
weights = .x$weights,
vB = 0.05,
dur = .x$frame_dur)))
Let’s see a few fits
map(modeldata[1:6,]$Logan_fit, plot)
[[1]]
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).
[[2]]
Warning: Removed 2 rows containing missing values or values outside the scale range (`geom_point()`).
[[3]]
Warning: Removed 2 rows containing missing values or values outside the scale range (`geom_point()`).
[[4]]
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).
[[5]]
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).
[[6]]
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).
Now let’s observe the resulting VT values.
Logan_outcomes <- modeldata %>%
select(sub, ses, Logan_fit) %>%
mutate(Vt = map_dbl(Logan_fit, c("par", "Vt"))) %>%
select(-Logan_fit)
ggplot(Logan_outcomes, aes(x=ses, y=Vt, colour=sub, group=sub)) +
geom_point(size=3, colour="black") +
geom_point(size=2.5) +
geom_line() +
scale_y_log10()
As we can see, they are reduced after blocking with calecoxib as expected.